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Abstract 

We introduce a new class of unfitted finite element methods with high order accurate numerical integration 
over curved surfaces and volumes which are only implicitly defined by level set functions. An unfitted 
finite element method which is suitable for the case of piecewise planar interfaces is combined with a 
parametric mapping of the underlying mesh resulting in an isoparametric unfitted finite element method. 
The parametric mapping is constructed in a way such that the quality of the piecewise planar interface 
reconstruction is significantly improved allowing for high order accurate computations of (unfitted) domain 
and surface integrals. This approach is new. We present the method, discuss implementational aspects and 
present numerical examples which demonstrate the quality and potential of this method. 

Keywords: numerical integration, level set, unfitted finite element method, isoparametric finite element 
method, high order methods, interface problems 


1. Introduction 

1.1. Motivation 

In the recent years unfitted finite element methods have drawn more and more attention. The possibility 
to handle complex geometries without the need for complex and time consuming mesh generation is very 
appealing. The methodology of unfitted finite element methods, i.e. methods which are able to cope with 
interfaces or boundaries which are not aligned with the grid, has been investigated for different problems. For 
boundary value problems with non-aligned boundaries methods such as penalty methods mm, the fictitious 
domain method Oil], the immersed boundary method [5], and other unfitted finite element methods [6] 
have been developed. For unfitted interface problems extended finite element methods (XFEM) have been 
developed in (among others) [71151 IHl [13 [H] ■ Also partial differential equations on surfaces which are not 
meshed have been considered using the trace finite element method (TraceFEM) [T^ . 

In the community of unfitted finite element methods mostly piecewise linear discretizations are con¬ 
sidered. One major issue in the design and realization of high order methods is the problem of numerical 
integration on domains which are only prescribed implicitly, for instance as the zero level of a scalar function, 
the so-called level set function. The use of standard integration rules (which ignore the cut position on cut 
elements) is not a good option as the integrand does not provide the necessary smoothness. 

The purpose of this paper is the presentation of a new approach which allows for high order numerical 
integration on domains prescribed by level set functions. The approach is new, robust and fairly simple 
to implement. The method is geometry-based and can be applied to unfitted interface or boundary value 
problems as well as to partial differential equations on surfaces. 

1.2. Literature 

We briefly discuss the state of the art in the literature to put the new approach into context. One 
approach that is often used in order to realize numerical integration on implicit domains consists out of two 
step: the approximation of the interface with a piecewise planar interface and a tesselation algorithm to 
divide the subdomains of a cut element into simple geometries, on which standard quadrature rules can be 
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applied. A famous example for quadrilaterals and cubes is the marching cube algorithm m- For simplices 
a detailed discussion of this approach can be found in (among others) [Ml Chapter 5], [15] for triangles and 
tetrahedra and in |161 HT] also for 4-prisms and pentatopes (4-simplices). Many simulation codes which 
apply unfitted finite element discretizations, e.g. [n muni mi m] make use of such strategies. In order to 
increase the accuracy of this integration one often combines this approach with subdivisions and adaptivity. 
Especially on octree-based meshes this can be done very efficiently [53]. However, this tesselation approach 
is only second order accurate (w.r.t. the finest subtriangulation) which complicates the realization of high 
order methods. 

One approach to solve this problem is the tailoring of quadrature points and weights which provide high 
order accurate integration rules for implicit domains. Such a construction of points and weights is based 
on the idea of fitting integral moments, cf. [M] [^. Although this results in accurate integration rules it 
has two shortcomings. First, the computation of quadrature rules is fairly involved. This aspect is typically 
outwayed by the resulting high accuracy. Secondly, the major problem, quadrature weights are in general 
not positive. This can lead to stability problems as positivity of mass or stiffness matrices in finite element 
formulations can no longer be guaranteed. 

For special cases satisfactory answers to the question of high order numerical integration strategies have 
been found which allow for an implementation of high order unfitted finite element methods. We mention 
a few interesting approaches. For quadrilaterals and hexahedra in |26j a numerical integration algorithm is 
presented which can achieve arbitrary high order accuracy and guarantee positivity of integration weights 
at the same time. The approach is based on the idea of interpreting the interface as a piecewise graph over 
a hyperplanes. In m an unfitted boundary value problem is considered. Instead of aiming at a high order 
accurate geometrical approximation of the boundary a correction in the imposition of boundary values is 
applied in order to recover a high order method. For the discretization of partial differential equations on 
surfaces, in [28] a parametric mapping of the interface from a piecewise planar representation to a high 
order representation based on quasi-normal fields is applied. Although the approach can not be carried over 
straight-forwardly to the case where also the integration on sub-domains is important, that approach has 
important similarities to the approach presented in this paper. 

In the literature of extended finite element methods (XFEM), there exist other approaches which are 
based on a tesselation approach and aim at improving the resulting piecewise planar approximations by 
means of a parametric mapping of the sub-trianguation [291130] . These approaches are typically technically 
involved, especially in more than two dimensions. Moreover, robustness of these approaches is not clear. 

The approach presented in this paper is similar to the approaches in [551HH130] in that it is also based 
on a piecewise planar interface which is significantly improved using a parametric mapping. The important 
difference is, that we consider a parametric mapping of the underlying mesh rather than the sub-triangulation 
or the interface. According to the mapping of the mesh the use of isoparametric finite element is natural. 

At this point, we would also like to mention the paper [3I| where the construction of a mesh deformation, 
which is also used in combination with isoparametric finite elements, is specifically designed to align a given 
mesh to a given interface position. The goal in that paper is to obtain a mesh that is conforming to a given 
interface without changing the mesh topology while keeping the quality of the mesh. Our approach is in a 
similar virtue. 

1.3. The concept 

In this paper we present a strategy for the numerical approximation of domains implicitly defined by an 
approximate signed distance function (f, the level set function. The approach is based on the assumption 
that a numerical strategy to deal with interfaces stemming from a piecewise linear level set function exactly 
is available. 

In the context of unfitted finite element methods integrals of the form fg f dx have to be computed for 
S being an interface or a domain which is only implicitly defined through the level set function (f. Consider 
the case S' = H n {^ < 0} for a polygonal domain H which has a suitable triangulation T. As an exact 
evaluation of the integral f dx is in general practically not feasible we consider the approximation of fli 
with the domain 

:=nn{ih(f <0} 
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where Ih is the continuous piecewise linear interpolation of 4>- To an explicit representation can be 
found on which quadrature rules can be applied element by element: 



fdxfv 




( 1 . 1 ) 


Here uji und Xi are the integration weights and points which depend on the cut configuration on T. The 
accuracy of this approach is limited by the approximation quality of fli^h which is only second order accurate. 
By a transformation of the mesh which is represented by an explicit mapping dt/j, parametrized by a 
finite element function, the piecewise planar interface is mapped approximately onto the zero level set 
of a corresponding high order accurate level set function, cf. Figure 0 



Figure 1.1: Main idea: The piecewise planar interface on the initial mesh is mapped close to the desired interface via a mesh 
transformation 'I';,. 


According to this mapping we have that is a high order accurate approximation to fli which, 

by construction, still has an explicit representation. The integral can hence be approximated as 





fdxKi EE w, |det(^'/,(a;*))| fi'^h{xi)), 

TgT i 


( 1 . 2 ) 


where the integration weights and points are the same as in (1.1). In contrast to (1.1) the accuracy of the 
quadrature in (1.2) is no longer bounded to second order accuracy but essentially depends on df/i. The finite 
element spaces (used in a discretization of an (unfitted) PDF problem) have to be adapted correspondingly 
which renders a resulting method an isoparametric finite element method: Let Wh be the finite element 
space corresponding to the piecewise planar interface approximation with Ti := {Ih4> = 0}. Then, after 
transformation with the appropriate isoparametric finite element space is 




Note that this implies that shape functions are not necessarily polynomials on the mapped domain. 

From a computational point of view only two things change compared to the case of a piecewise planar 
interface. First, a suitable (approximated) mapping has to be constructed. A major contribution of this 
paper is the discussion of this. Secondly, the deformation has to be considered in the definition of the shape 
functions and the quadrature of deformed elements. The treatment of the latter aspect is well-known and 
is not discussed here. 

Owing to the combination of a tesselation algorithm for a piecewise planar interface and a parametric 
mapping, we obtain the four most important features of this approach: 

• An explicit high order accurate geometrical approximation of the exact interface. 

• Guaranteed positiveness of quadrature weights on the interface and in the sub-domains. 
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• Unchanged cut topology compared to the piecewise planar case. 

• Easy integration into existing codes as the approach builds on a piecewise planar interface. 

Together, this renders the method highly accurate, robust, efficient and fairly simple to implement. We 
again note, that the improved geometry approximation is based on a parametric mapping of the initial 
mesh, s.t. a combination with isoparametric finite elements is crucial. 

1 . 4 . Content and structure of the paper 

The main purpose of this study is the presentation of the new approach which allows for an explicite 
high order geometry approximation of domains which are implicitly defined as level sets. We focus on the 
discussion of an efficient construction of a suitable mesh deformation and related implementational aspects. 
Although this discussion goes beyond simple concepts in that it addresses important details of the method, 
we do not aim at a thorough error analysis, yet. This is topic of ongoing research and will be published in 
a different paper. 

The paper is structured in the following way. The starting point for the method is the ability to 
deal with interfaces which are piecewise planar. As we mainly consider simplicial meshes this piecewise 
planar interface typically corresponds to a (continuous) piecewise linear level set function. We consider a 
corresponding decomposition strategy as given and do not discuss it here. Instead, we refer to the literature 
instead, cf. [HI Chapter 5] or [T3 Chapter 4]. 

The crucial component of the new approach is the construction of an explicit mapping 'P/j which is suitable 
for implementation. The construction of '^h is discussed in section and consists of several steps. We first 
characterize the desired properties of the transformation ’^h for a general case. Afterwards we restrict to 
an important special case with simplicial elements and a given level set function which is an approximate 
signed distance function. For this case we present an explicit construction of a suitable mapping. 

Numerical examples demonstrating the quality of the explicit geometrical representation obtained by 
this new approach are presented in section and section While the examples in section focus on the 
accuracy of the geometry approximation, in section a high order unfitted (isoparametric) finite element 
formulation for an interface problem is considered which proves the practical use of the method. 

2. Construction of the mesh transformation "ifu 

We start with the formulation of properties that we demand from a suitable transformation Let 
Ti be a piecewise planar interface and fLe space of continuous vector-valued piecewise polynomial 

functions of degree k. We seek for a transformation G \VhY-i s.t. 

4'^ = argmin dist(T/j(ri),r) or at least dist(4'/j(ri),T) < (2.1) 

'I'hGV 

with V* C fLe subset of admissible transformations which fulfil the following constraints: 

• Homeomorphy: 

e [C°(U)]‘^ (2.2a) 

• Shape regularity: Consider an arbitrary element T in the triangulation T. Let be the transformation 
from the reference element T to T = $/j(T). Then the mapped element T = '^h(T) = (dt/j o $ft)(T’) 
should also be shape regular. Moreover, the shape regularity should be comparable to the shape 
regularity of T. We therefore define the transformation 


and demand 

«:(V0,,) < ^^(V®,,). (2.2b) 


for a given constant C > 1. 
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Locality: 


'^hix) ^ X only in the vicinity of the interface Li. 


(2.2c) 


The last constraint is desirable for efficiency reasons but not crucial. 

Note, that due to the fact that Ti is already a second order accurate approximation to the interface, the 
transformation is essentially only a local high order correction. This is in contrast to the transformations 
tailored in j3lj . Further note, that 'I'/i is obviously not uniquely defined with the requirements formulated 
in and ( |2.2a[ )- p.2c| ). Different choices are possible. For the case of a simplicial mesh we introduce a 
special choice for 'I'ft, in the following. 

In order to construct a suitable transformation ^'ft the above given characterization is too general to 
be of direct practical use. In the remainder of this study we restrict to simplicial meshes and assume the 
interface to be described by an approximate signed distance function (j) which is a continuous piecewise 
polynomial function. Further, the piecewise planar approximation Fi is assumed to be the zero level of the 
piecewise linear interpolant Ih4>- These restrictions allow for an explicit construction of a function d^ft which 
is of practical use. 

we 


We derive this in several steps. In section pTl we introduce notation and assumptions. In section 2.2 


give an ideal mapping 'I' which fulfils 'l'(ri) = Ffc := {(j) = 0} (and hence (2.1)) and is defined pointwise. 
Together with a suitable modification (section 2.3) and an approximation of this function in (section 

2.4) we implement conditions (2.1), (2.2a) and (2.2c). In cases where the interface is well resolved by 


the piecewise linear level set function Ih(t> the transformation 4('ft is close to the identity and the shape 
regularity of the transformed mesh is ensured. Nevertheless, for practical applications a mechanism to 
ensure robustness even if the interface is not well resolved becomes necessary. The aspect of controlling the 


shape regularity as in (2.2b) is implemented with a limitation step which is presented and discussed in detail 
in section 12.61 

In section [2.7| we summarize the algorithm to determine Tft and the most important properties of the 
method from a computational point of view. 


2.1. Notation and assumptions 

We introduce some basic notation and assumptions. D is a polygonal domain in d = 2, 3. It is 
decomposed into a shape regular partition T of D consisting of elements T which are simplices. Inside Q an 
internal interface F (F n 912 = 0) separates two domains and ^ 2 . The interface and the subdomains F, 
rii , f 22 are implicitely described by a level set function with F = {a:| = 0}, fli = {a:| < 0} 

and ^2 = {x\ (j)’^^{x) > 0}. The level set function is a signed distance function: 




dist(a:,r), a; S fli, 
—dist(a:,r), xGfl2. 


We make the following assumptions on the smoothness of the interface F and the resolution of F by the 
triangulation T. 

Assumption 1. We assume that the interface F is -smooth, m > 2. Then there exists a S > 0, such 
that (j)^^ G with := {a;| \\x — y \\2 < S for a y GT}. 


Assumption 2. Let h he the maximum mesh size of the triangulation T. Then we assume that there holds 
h < S with d as in assumption^ 


Assumption]^ can always be satisfied for sufficiently small mesh sizes. Let denote the space of continuous 
piecewise polynomials of degree k. The level set function is approximated by a function (j) G V^, 
1 < k < m — 1. We think of ^ as a projection of under a suitable projection operator : L^{Ll) -G Vjf. 
In practice one often only knows (j) and not We add an assumptions on the approximation error induced 
by the projection 11'^. 


Assumption 3. With (j) = and Ffc := {(j) = 0} we assume that there holds 
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(2.3a) 

(2.3c) 


dist(rfe,r) < ch'^+^. 


(2.3b) 

(2.3d) 


l|V(/>||oo.o^ <c||vriloo,n^, 

with constants c only depending on k and F. 

Note that assumption is justified for IF^ the standard projection on V, or the projection discussed in 
section 12.51 

The standard nodal interpolation of (p into the space of piecewise linears is denoted by Ih4>- The zero 
level of that level set function Fi := {Ih<P = 0} is piecewise planar which allows for an explicit representation 
of the interface and the sub-domains. Fi is a second order accurate approximation to F. Let 

T’" := {T e T| T n Fi 7 ^ 0} and :={x£T\T e T^} 

be the set and region of cut elements, respectively. Due to assumption the interface is resolved in the 
sense that D'" C With the assumed smoothness (m > 2 in assumption^ there also holds: 

\\Ih(t> - (pWL-^iOP) < c/l^||())®"'|| 772 (nr). (2.4) 


We further introduce the following notation for the extension of cut elements by its neighbors: 

Tr.+ := {T e T| 9r n 7 ^ 0} and := {x S T] T G T^’+}. 


2.2. Ideal transformation 4' 

In this section we define an ideal transformation 'F which is defined pointwise and is not a finite element 
function. We assume the ideal case that (p = and that the assumptionsandare valid. This specifically 
implies that Wp is continuous in the set of cut elements D'" and F = F^. At the end of the section we discuss 
the case p ^ p^^. 

We ask for 4t(Fi) = F^ = F. This only determines the transf orma tion tk on Fi. Hence, we use a condition 
which describes a reasonable extension on D'". Later, in section 2.4 the continuous transition to 4t;i(a;) = x 


away from the interface is implemented within the approximation in a finite element space. 


Definition of the mapping We choose a transformation dt which maps iso levels of the piecewise linear 
approximation of the level set function {Ihp = c} onto the corresponding iso levels of the level set function 
{p = c}, i.e. we ask for 

Ihp = po^ V a; G (2.5) 


Such a transformation is sketched in Figure 2.1 Note that (2.5) implies 'I'(Fi) = F. The mapping can also 



Figure 2.1: Ideal transformation: Piecewise linear level sets are mapped on exact level sets. 


be formulated as the following problem: For every point x € ^2^ with (linearly interpolated) level set value 
c = Ihp{x), find a point x* in D such that p{x*) = c. 

Note that (2.51 does not define dt uniquely, yet. In order to determine a unique transformation, for every 
point X in D'" we search for the corresponding point x* G {p{x*) = c} with smallest minimal distance to x. 
Due to a: G and assumption this point is unique. With the help of the signed distance function p we 
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can characterize x* = x + rs using the search direction s = V4> and the distance r. By specifying the search 
direction s{x) for each point in we can define dt as 

•.= x + r{x) s{x), x G with r(x) G K such that (2.6a) 

Ih(l>{x) = 4>{x + r{x) s{x)). (2.6b) 

With (j) = we have that s = is continuous in which results in € [(7°(r2^)]‘^. By construction 

we have 'f'(a:v) = xy for every vertex xy in the mesh as there holds Ih4>{xy) = 4>(xy). 


The non-ideal case (j) ^ In practice one mostly has cj) ^ cjf^. Then, typically s = V(() is not continuous 

any more and as a result the same holds for a correspondingly defined dt. There are essentially two ways 
to deal with this problem. Either continuity of the search directions is restored by means of a projection 
n® of W(j) into Possible projections are for instance an L^(r2^)-projection or the projection discussed 

in section [2(5| Alternatively, one could stick with the discontinous search direction s = Wcj) accepting the 
resulting discontinuous This is acceptable as finally, for a realization of the transformation, 'k is projected 
into a continuous finite element space with a projector B’^ similar or identical to 11®. In our experience both 
approaches give similar results. Nevertheless, in the remainder of the paper we consider the use s = n®V(/) 
as it simplifies the discussion of the method. We also mention that it is not necessary to have a high order 
accurate approximation of for the search direction. In the numerical examples in sect ion [3| and section 
l^a continuous piecewise linear search direction (using the projection discussed in section 2.5) is used. 

Another important aspect of the non-ideal case (f) ^ is the fact, that the transformation is typically 
not smooth (only Lipschitz continuous) within each element T as 4> typically has discontinuities in all 
derivatives across element interfaces. A typical situation leading to a non-smooth dt is sketched in Figure 


2.2 In section 2.3 we propose a modification which (alongside its other benefits) results in a transformation 


which is smooth within all elements. 



Figure 2.2; Sketch of a situation where the ideal transformation is not smooth within a given element (^ 2 ). Points on the 
piecewise planar level set (green) are mapped on corresponding points on the corresponding level set of </< (blue), if inherits 
the missing regularity of (p with a shift in locations: kinks across element interfaces in p result in kinks of ili within elements. 


Evaluation o/tk.'. At each point where the transformation tk has to be evaluated, the solution of problem 
(2.6b) for a fixed x has to be implemented. Let T be the simplex, such that x G T. The search direction s(a;) 
and the value (Ih4’){x) can easily be evaluated at the given point. This information is available element- 
local. In order to find the corresponding projected point x*, it seems natural to try the following Newton 
iteration: 

(4(/>)(a:) - <(>(x'=) 


X = x, 




= x^- 


V(^{X^) * 5(x) 


s{x) 


(2.7) 


This procedure has two disadvantages. First, (j) may have kinks at element interfaces (as it is only piecewise 
polynomial) which may complicate the solution of the non-linear problem. This problem can be tackled 
by replacing the Newton-iteration by more sophisticated methods. More importantly, for € T it can 
happen that x^ ^ T, such that evaluations of (j) need to be carried out in a neighborhood of element T. 
The evaluation of finite element functions on neighbor elements is often troublesome, esspecially if one is 
concerned with an efficient parallelization. We propose a modification which solves both problems in the 
next section. 
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2.3. Localized (piecewise smooth) transformation 'I't 

We modify the previously discussed construction of the deformation in order to render an implementation 
less complicated and more efficient. We achieve this with a localization of the computations on elements. 
That means that we avoid evaluating finite element functions from neighbor elements. To this end we first 


introduce a discontinuous (across element interfaces) transformation 'I't in section 2.3 which is afterwards 
projected into the space of continuous finite element functions with a projection such as the one discussed 
in section 12.51 

We define a modified transformation by its element-wise contributions and only consider elements 
T . On each element T the level set (f is a polynomial of order k. Let £t be the extension operator so 
that £t{4>) denotes the polynomial on the restriction of which on T coincides again with 4>. As (f is an 
accurate approximation of a smooth function, is a good approximation to (f also in the neighborhood 


of T. Replacing (j) with £t{4') in (2.6b), results in the definition of 'I't: 


'f' 7 ’(a;) := X + r^ix) s{x), x GT with rx{x) G 
Ihfiix) = {£t{(I))){x + rrix) s{x)). 


such that 


As the function to evaluate in (2.8b) is now a polynomial on the Newton iteration 

{Ihffix) - {£T{(t>)){x^) 


^k+i = _ 


V(fT(</>))(a;'=) • s{x) 


s{x) 


(2.8a) 

(2.8b) 


(2.9) 


is much simpler to evaluate and should converge even faster. For moderate order k even exact root finding 
algorithms could be applied. Another advantage of this approach is that we can guarantee that the resulting 
transformation is smooth within each element, cf. Figure [2^ 



Figure 2.3: Sketch of the same configuration as in Figure r2.2| this time for the localized transformation 'I'j’. Points on the 
piecewise planar level set (green) are mapped on corresponding points on the corresponding locally extended (dotted) level set 
of b (blue) resulting in a smooth function iky within each element T. 


We briefly explain why we expect this to be a valid approximation of the original problem. To this end, 
we assume that the interface is resolved. In that case £Ti4‘) only has to be evaluated in a small neighborhood 
of T, denoted as oj{T). In that neighborhood a simple triangle inequality gives 

\\£t{()) — </'||oo.w(T) < — <(>®^||oo,w(T) + Wf’ — ()‘^^\\oo,u:{T) (2-10) 

where the latter term can be bounded in due to assumptionand the first term on the right hand 

side can also be bounded in (!I(/i^+^) with a Bramble-Hilbert argument. 

'I't is discontinuous across element boundaries. However, the jumps over element boundaries should be 
in the order of Again by a projection H'^ we make sure that the final discrete transformation 


is a continuous finite element function. A simple and efficient choice for H'^ is discussed in section 


2.5 


2 . 4 . Approximation with finite elements 

Until now we have discussed how to define and evaluate 'h in Our goal, however, is the construction 
of a finite element function 'L/j G which fulfils condition (2.1), but furthermore guarantees (2.2a) and 


















We write '^h G [V^]'^ as = a; 


dfi with dh G WhY function describing the deformation of the 


mesh. Now we can split 


iVp 




dh = di 




i=l 


i-1 


Vj¥’?(x) 


with Ui G the degrees of freedom corresponding to the basis functions (p^ which are located at the 
interface, supp((/5f )nn^ ^ 0, and Vj the remaining degrees of freedom (supp(i^^) = 0). By construction 
we then have d'/jjQr = (a; + d)i)|nr and d^|nr = 0, i.e. d^ has no contribution in 

This splitting decomposes the problem again into two parts. First, has to be chosen such that '^h is a 
good approximation to da in Secondly, has to be chosen such that = a; in 0 \n^’+. The second 

problem immediately suggests to choose d^ = 0. Hence, (2.2a I and (2.2c) are fulfilled and we only have to 
consider the approximation in (2.1). On that account we apply a projection H'*' of (or da) into 


This projection could be a nodal interpolation (only if da is used), an L^(f2'")-projection or the projection 
in section [2.5| and we expect the following assumption to be valid. 


Assumption 4. Consider da* € {da, da^}. With da/j = H'^d'* there holds 


1^/* - < ch'=+l||da*||^.+i(Tr) < ch’^+^l^Y 


( 2 . 11 ) 


where the last inequality holds due to da*(x) x. Further there holds 

ll'E'a.-da*||^or <ch'=+i. 


( 2 . 12 ) 


A rigorous justihcation of this assumption (for both, da and daf) requires further analysis which is topic of 
ongoing research. 

Consider da^ = H'^d'. Dehning T/i := da^(ri) we expect the following error estimate to hold: 

dist(r;„r) < ll.^'^'lloo.r. = = \\r^o-^h-h(l^\\oo,Ti 

< ||(^®’'od',i-())od',j||oo.ri+||0od'-())od',j||oo,ri + ||<?i’od'- 4(^||oo.ri (2-13) 

< U‘'^-^\Unr + |!V,^||^,or ||da-da,||^,r,_ <ch"+^ =° 


where in the last step we used assumptions and For the case da/j = H'^d'T the estimate still applies 
(replace da with dai) except for the term ||(/)odaT — Ih4’\\oo,ri in (2.13) which is not exactly zero but can be 
estimated as follows: 


||()) odaT — //t(/>||oo.ri < max { \\£T[(t)) odai — -f/i'ii’||oo,rinT +\\^t{Y) odaT||oo,T} 

TGTr s_„_/ 


=0 


< -'P\\oo,ui(T) ch 


fc+i 


2.5. An Oswald-type projection 

At several occasions throughout this paper we mentioned a generic projection operator H : [L^(H*)]"’ —>■ 
[Vjf]" with H* = H or H* = H'" and n € (1, d}, for instance H"^, H® and H'*'. An obvious choice for all these 
projection is a simple L^(H*) projection. However, to avoid the solution of global linear systems it is worth 
mentioning a simple and more efficient alternative which we also used in the numerical examples. For ease 
of presentation we only consider the case where we seek a projection H : L^(r2'") —>■ V^, i.e. the scalar case 
on the interface region. The translation to H and/or vector-valued functions is then obvious. 
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The Oswald-type interpolation operator consists of two steps: A projection into a discontinuous finite 
element space and an averaging into Vjf. Let := {u\t G 7^^(T)| T S T^} be the space of piecewise 

polynomials which are discontinuous across element boundaries. Due to the missing continuity-restriction, 
the standard L^-projection : L^(n) —>■ can be computed in an element-by-element fashion and 

can be implemented very efficiently. Next, we define a simple averaging operator IL^'' : which 

is a high order version of the Oswald interpolation operator [32]- Let Wi := {T G T’"! supp((/ 7 i) n T 7 ^ 0} 
be the set of elements where the ith basis function ipi of Vjf is supported in Further, let Ui^r be the 
coefficients such that 

u = ^^\T, u G 

i Teui 

We then define 11 ®-''(m) := ^ - Wiipi 


with Wi = —— Ui^T if supp((/ 3 i) n o'" 7 ^ 0 and Wi = 0 otherwise. 

Teuji 


Together, we obtain the projection operator 11° : L^(0) —>■ Vjf with 11° = 11®'' o We note that in an 

implementation of 11 ° the finite element space V'^’‘'*®° does not need to be constructed explicitly. 


2.6. Shape regularity 

So far we assumed that the interface is smooth and well resolved. An approach which aims to be of 
practical use should however be capable of dealing with non-smooth or underresolved interfaces. In this 
section we characterize a simple sufficient condition to ensure shape regularity. Based on this we discuss two 
different cases. First, we assume the nice case, where the interface is well-resolved and smooth and show 
that shape regularity is not an issue. Secondly, we discuss the case where the interface is not resolved and 
propose a simple strategy to ensure shape regularity. 

Before we discuss the details of both cases, we introduce transformations between the reference simplex, 
the simplex before and after the transformation with 'I'/j and an intermediate curved simplex. 


Simplex transformations:. A planar simplex can be characterized as a reference simplex under an affine 
transformation. The mesh obtained after the transformation with can be interpreted as a collection 
of simplices under the concatenation of two transformations: an affine one and another one parametrized 
by a finite element function. This characterization facilitates the investigation of the shape regularity of 
the transformed elements. We therefore introduce notation for the (curved and planar) simplices and the 
transformations relating them. In figure 2.4 the simplices and transformations are also depicted. 

By T and x we denote the reference simplex and its coordinates. The affine mapping transforms the 
reference simplex to the planar simplex T = $^(T). 


X = <i)(a:) = A - X xq (2.14a) 

The final (curved) element T is obtained after the finite element transformation of T by The shape 
functions (corresponding to 'I'/i) are typically defined w.r.t. to the reference simplex T, s.t. we have 

y = 'f>h{x) = X + '^di{ipiO ^f^){x) xo + A-x + '^di- ipi{x). (2.14b) 

i i 

Thus, we have T = &h{T) with 0/i = 4'^ o We can also characterize Qh as a mapping which first 
applies the non-planar deformation and the affine mapping afterwards. We therefore introduce the mapping 

s.t. 

0;, = (4',, o $;,) = ($^ o 4-^) (2.14c) 

with 

y = i/h{x) = x + '^di- ipi{x), dt=A~'^-di. (2.14d) 
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Figure 2.4: Reference and transformed domains before and after curving and before and after affine transformation to the 
physical domain. 


Let xv and xy be the coordinates of the vertices of T and T. Then 'I>ft,(xy) = xy and d>ft,(xy) = xy. 

From the previous observations and we can conclude that 

V0 = V<^h = A. (2.15) 

i 


Hence, if the transformation ■ T ^ T (curving of the reference element) and the transformation : 
T ^ T (affine transformation) are well-behaved, the same also follows for 0 : T —>■ T. 

As a measure to quantify the shape regularity we consider the relative spectral condition number of the 
Jacobian. We have the simple estimate 

k(V0) < k(V$,,) • (2.16) 


We may assume that the initial mesh is shape regular such that K(V$ft) < C for reasonable moderate 
number C. In order to fulfil (2.2b) it remains to ask for the impact of the non-linear transformation, i.e. 
K(VTft). 

We have = / -I- Wdh with dh{x) the deformation of T w.r.t. the reference simplex T. As dh is a 

polynomial on the reference triangle, we can bound its gradient 




oo^T 


(2.17) 


with llalloo f = max^gy 


g(x )||2 for a vector-valued function q and HQIlg^ f = max^g^ ||(5(x)||f for a matrix¬ 
valued function Q. The constant Ck only depends on the polynomial degree. Due to the choice of norms, it 
is sufficient to have (2. for a scalar-valued polynomial. 

A simple sufficient condition for the shape regularity, formulated as K(V'I'ft^) < C for a chosen constant 
C > 1, can be derived as follows. With elementary calculations we have 


l | Vd /^||2 

ll(V'i'^)-'ll2 


< i + I|V4||2 
<(1-||V4||2)-i 


< 1 + ||V4||;^ 
<(l-||V4||f)-i 


<1 + I|V4L 

<(i-l|v4lL,t) 


-1 
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With (2.17) we can easily conclude that 


is a sufficient condition for 


\\dh\\ 


oo.f — 


C-1 

Ck{C + 1 ) 


< C. 


(2.18) 


The resolved case:. Due to \\Ih4' ~ 4>\\L'^(nr) < ch? there holds dist(ri,rfe) < ch? and we can conclude that 
lld/illj^ f < ch. This implies that the deformation gets arbitrary small which implies (2.18) for sufficiently 
small h. We can thus conclude, that the shape regularity is easily provided for sufficiently small mesh sizes, 
i.e. for a sufficiently fine resolution of the interface by the mesh. 


The under-resolved case:. The considered construction of the transformation may lead to undesirable de¬ 
formation, e.g. self-intersecting geometries and arbitrary small angles if the interface is not sufficiently well 
resolve. This may lead to a blow up in K(\7^>h) or even a singular at some points. 

To avoid these situations, we introduce a “barrier step” which garantuees the quality of the resulting 
mesh. The idea is simple. We replace the deformation d = — x (or d = 'I't — x) with a limited deformation 

d = minly/i • s/||s||, d}. (2-19) 

On a fixed element T, we then have d^ = (V<i>/i)“^n'^d and hence 

||4||oo < ctC^7 (2-20) 


a constant depending on 


with Ctt the L°° stability constant of the projection 11'*' and ct = h||(V$ft,)“*| 
the shape regularity of T and quasi-uniformity of T*". 

For sufficiently small parameter 7 we can then ensure condition (2.18) and get shape regularity indepen¬ 
dent of the resolution of the interface. The price, however, is a reduction to an essentially only second order 
accurate approximation if the limitation step has to be used {d d). 

In cases where the interface is resolved the limitation is not needed which leads to d = d. In the 
worst case, where the interface is not resolve, one essentially recovers the quality of the piecewise linear 
approximation. In both cases the deformed elements are shape regular which results in a robust method. 


2. 7. Summary of computational aspects 

To clarify on the (low) computational complexity of the resulting scheme we briefly summarize the 
algorit hmic structure for the computation of We restrict to the modified version of the algorithm, cf. 
section 2.3 which is also used in the numerical examples. To this end we determine the coefficients di G 
of the representation 'I';i(x) = x -f di(pi{x), di G where ipt are the basis functions of V)f: 


1. Set coefficients and counter (for later averaging) to zero: 
di = 0 and = 0 , i = 1 ,.., dim(V)f). 

2. Compute element contributions ( 11 '*'®''d/^): 

Loop over elements T gT^: 

(a) Determine the 'I'/i|t G [7^^(T)]'* by solution of the local L'^iT) problem: 


'^h\T-'vdx= / 'I'T-'cdx 'ivG[P^{T)\‘ 


For each integration point Xj this involves: 
• The evaluation of 


This requires the solution of (2.8b) by the Newton iteration (2.9). 

• The limitation step: 

If the deformation 4'r(xj) —x^ is too large, it is restricted according to (2.19). 
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(b) Add element contribution to global coefficient vectors (for all d.o.f. i of T): 

• di di + df where df are the basis coefficients corresponding to d'/ilr 

• Tii m + l 

3. Average the element contributions ( 11 ®'''): 
di di/rii if > 0 

In this overview we neglected the computation of the search direction s = n®V0. The corresponding 
algorithm goes the same lines with the only difference that the evalution of the right hand side functional 
in the (element-)local projection is much simpler. 

The resulting deformation 4'^ has the following properties: 

• Only elements in the neighborhood of T are deformed. In \ 17'"’+ the mesh stays unchanged. 

• Where the resolution of the interface is sufficiently high, the interface can be resolved with high order 
accuracy. 

• In any case shape regularity can be guaranteed. 


3. Numerical examples I: Geometry approximation 


Preliminaries 

In the following we consider complex geometries described by level set functions and their approximation 
with the approach proposed before. For the barrier parameter we choose 7 = 0.1 (independent of the degree 
k). To solve the one-dimensional nonlinear problem in (2.6b I we choose a (relative) tolerance of 10“^'* which 
resulted in only of 2 — 3 iterations for each point. 

For the quasi-normal held s we use a continuous piecewise linear held. In all examples we use the 


transformation = n°('I'T), cf. section 2.3 The implementation is included in an add-on library ngsxfem 
to the hnite element library NGSolve 

We are interested in the geometrical error between the discrete interface F/j = 4';i(ri) (with an explicit 
description) and the ideal interface F = {(j)‘^^{x) = 0} (with an implicit description). Note that the error 
consists of two approximations: The approximation of the level set function (j) « and the approximation 
of Ffe = {(() = 0} with 4';i(Fi). The error we are interested in is the maximum distance between the discrete 
and the exact interface, dist(F;i,F). 


3.1. A two-dimensional example: The flower-shape 

The geometry in this example is inspired by the one in [53] . Due to its shape we call it the “hower-shape” 
example. The domain is a square domain 17 = [—1,1]^ and the interface is given as 

F := {(a;i, 0 : 2)1 xi = R{0) cos{0), X 2 = R{0) sin(0), 9 e [0, 27r]} (3.1) 

with 

R(9) = rg -h 0.1 sin(a;0), rg = 0.5, uj = 8. 

F is the zero level of the level set function (j){x) = ^/xl -\- X 2 — R{9) with 9{x) = aictaYi{xi/X 2 ). The inner 
domain is Di := {x € D, (j){x) < 0}. A sketch of the geometrical configuration is displayed in Figure [0] (left) 
together with the initial unstructured mesh which consists of 230 triangles. Starting from that mesh, local 
rehnements around the interface are successively performed 12 times. The final mesh consists of roughly 
1.5 million triangles. We consider different polynomial degrees k, where we always use the same polynomial 
degree for the approximation of (j) and the deformation 4'/i. 

In this example f is not a signed distance function, s.t. dist(Fft,,F) f ||<(>||oo,rfe- However, there holds 
||V0|| > 1, s.t. we have dist(F;i,F) < In Figure 

is an upper bound for the distance. 

If the mesh size h in the vicinity of the interface is sufficiently small the geometrical error decreases 
with the (optimal) order fc -|- 1. For coarse grids, where the resolution of the interface is not sufficient the 
transformation has to be limited in order to garantuee shape regularity. On the corresponding grids we can 
not expect to observe a high order convergence which we also observe on the levels 0 — 2. 


3.1 


the convergence of 


is depicted which 
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Figure 3.1: Mesh and convergence behavior of geometrical error for the flower example. 


3.2. A three-dimensional example: The gyroid 

In this example we consider the cube = [—1,1]^ with an interface prescribed by the zero level of the 
following level set function: 

(j){x, y, z) = cos(7ra;i) sin(7ra:2) + cos(7ra;2) sin(7ra;3) + cos(7ra:3) sin(7ra:i) (3-2) 


A sketch of the interface and the initial unstructured mesh, consisting of 577 tetrahedra, is depicted in 
Figure 3.2 In a neighborhood of the interface F there holds ||V(^|| > 1, s.t. we again have that the 
computed quantity H^^Hoo.rh is an upper bound for dist(r;i,r). We consider seven adaptive refinements 
only towards the interface. The final mesh consists of around 40 million tetrahedra. Note that this time 
the interface is not contained in such that ^ fl. However, for the approximation of F and the 

evaluation of the error this is not important as (p is trivially extended to 



Figure 3.2: Mesh and convergence behavior of geometrical error for the gyroid example. 


The results in Figure reveal that the asymptotic behavior of the method is as good as in the previous 
two-dimensional example, the order of convergence is clearly visible. However, the pre-asymptotic 
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regime is significantly larger. This is due to the more difficult geometry. We observe that until refinement 
level 4 the error drops only with 0{h?). However, once the interface is resolved by the piecewise linear 
interface Ti, the limitation of the deformation is no longer necessary and the high order resolution of the 
interface kicks in. On the finest level, the accuracy between different polynomial degrees differs dramatically. 
An increase of the polynomial degree by one order results in 2 orders of magnitue decrease in the error. 

3.3. Summary of numerical examples in this section 

These example represent configurations with realistic and challenging complexity. We conclude that the 
method is highly accurate if the interface is smooth and well resolved. If the interface is not well resolved, 
we at least obtain a robust method with the “standard” second order accuracy. 


4. Numerical examples II: A high order unfitted isoparametric finite element method for an 
interface problem 

J^.l. Elliptic interface problem: The disk 

Problem description:. We consider the problem from m section 2.5.1.4], where an elliptic interface problem 
of the following form is considered. 


div (aVu) = / 

in Hi, 

* = 1,2, 

(4.1a) 

[aVu • n] = 0 

on F, 


(4.1b) 

II 

0 

on F, 


(4.1c) 

u=9D 

on 9H. 


(4.1d) 


Here Pi, i = l,2 are domain-wise constants, ( 01 , 02 ) = (2,1), {^ 1 ,^ 2 ) = (1,3/2). We consider a circular 
interface T with radius R = 0.6 as interface T in the square domain [—1,1]^. The data go and / is chosen 
such that the solution to (4.1) is 


u{x) = 


a2U{r{x)) 

aiU{r{x)) 


■ P2, X e Oi 
'Pi, X G fl2 


. nr 


with U{r) = “ 




The Nitsche-XFEM discretization:. We consider an isoparametric version of the discretization in [ni section 
2], where a combination of an extended finite element (XFE) space combined with a Nitsche formulation 
for the unhtted interface is applied. The discretization is compactly described next. For a more thorough 
discussion of the underlying (planar) discretization we refer to [iTl chapter 2]. 

Let TZi be the restriction operator to domain with fli^h ■= {IhP < 0} and fl 2 ,h '■= {Ih'P > 0}, such 
that (7?,jit)|n = Sij ■ u\q. with Sij the Kronecker delta symbol. The extended finite element space is then 
defined as © TZ 2 V^. Basis functions from are continuous within the sub-domain but 

may be discontinuous across Fi. According to the transformation we define the isoparametric finite 
element space = {vh o '^f^\vh G V)/} with correspondingly mapped shape functions. 

The finite element space is not conforming w.r.t. the interface condition (4.1c). To incorporate the 
interface condition we consider a Nitsche-type discretization: 


Find Uh G Vf such that 


a{uh,Vh) + Nh{uh,Vh) = f{vh) for all Vh G 


(4.2) 


with Nh{u, v) := Nf^{u, v) + Nf{v, u) + Nf{u, v) and 


a{u,v) := / aiVuVvdx, f{v) := / fvdx 

2—1 2 


Nh{u,v)'= l-aXu-n}lPvjds, 

dr. 


^h{u,v):=J^ ‘^^^IPujlPvjds, 
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with F/i = a = ^{ai+a 2 ) and A a stabilization parameter that only depends on the shape regularity 

of the mesh. 

The terms a{uh, Vh) and f{vh) ensure consistency inside the subdomains, Nf^{uh, Vh) ensures consistency 
w.r.t. the interface condition (4. N^(vh,Uh) render the discrete l.h.s. operator symmetric while keeping 
consistency and Ns{uh,Vh) ensures stability in a consistent way. 

An important detail in this discretization is the choice of the element-wise defined averaging operators 
We define the averaging weights Ki, s.t. ■= '^i'a|ni + K 2 w|n 2 - We consider the “Heaviside” choice 
where = 1 if \Ti\ > ^\T\ and Ki = 0 if |Ti| < ^\T\. Here [Til refers to the cut configuration on the 


undeformed mesh. This choice in the averaging renders the scheme in (4.2) stable for arbitrary polynomial 


degrees k, independent of the cut position of T. The interested reader is referred to the error analysis in [m 
section 2.3], specifically remark 2.3.1 therein. 

This test case has been considered for k = 1,2,3 in na section 2.5.1.4] with a piecewise planar approxi¬ 
mation of the interface (on a eight times adaptively refined subtriangulation). We now consider k = 1, 2,3,4 
in V[ in combination with a level set function and a transformation of the same polynomial degree. 

Numerical setup:. We consider the mesh displayed in Figure o as our initial mesh and successively apply 
uniform refinements of the mesh and measure the error in the norm, the broken (semi-) norm as well 
as the error in the interface condition (4.1c) in the L^(r/j)-norm. To solve the arising linear systems we use 


a sparse direct solver. We note that the arising linear systems are extremely ill-conditioned. We observed 
that in some cases, for A: > 2, the linear systems were so ill-conditioned that the initial residual of the linear 
system could only be reduced by a factor of 10“®. 

Remark 4.1. Efficient linear solvers for this kind of problem are rarely addressed in the literature, especially 
for high order unfitted finite elements. The only robust strategy (without additional stabilization) for linear 
systems arising from Nitsche-XFEM discretizations, that we know of, which is based on a rigorous analysis, 
is presented in m- In that paper however only the case k = 1 is considered. The design and analysis 
of linear solvers for linear systems arising from Nitsche-XFEM discretizations remains a challenging task 
which requires further attention. 


Numerical results:. We note that no barrier step (cf. section 2.6) has been necessary after one refinement. 
Hence, we expect to be in the asymptotic regime w.r.t. the geometry error such that U^^Uoo.rh < 
holds. We observe the (optimal) convergence rates for the volume error that are predicted in the error analysis 
in [HI section 2.3]. Note that error sources due to geometry approximation and numerical integration have 
not been considered in that error analysis. The error in the interface condition |/3u] = 0 converges with 
0(/i^+^) which is even half an order better than predicted. 


“-“hllryn) 11“ — “hll/ryoiuna) ll[/3“/illUyrh) 



refinements (uniform) refinements (uniform) refinements (uniform) 

Figure 4.1: Error convergence in different norms for the circle example. 


16 















Conclusion 


A new unfitted finite element method with a high order geometrical approximation of level set domains 
is presented and discussed in detail. The method is efficient and easy to implement. Numerical examples 
reveal its potential in handling complex geometries robust and highly accurate. 

As an example for the potential of the method, the method has been combined with a rather standard 
Nitsche-XFEM discretization for an unfitted interface problem and optimal order of convergence has been 
observed. The potential of the method goes beyond this specific example. As the method is geometry-based 
it allows to consider high order methods for a large range of unfitted finite element methods. However, 
many open questions remain. Stable discretizations and suitable linear solvers for high order unfitted finite 
element methods are difficult to develop and analyze. 

The presented method requires further attention. A rigorous error analysis of the presented method is 
missing, yet. We will address this in a forthcoming paper. Further, extensions to non-simplex meshes and 
time-dependent problems (with moving interfaces) are interesting topics for future studies. 
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